SUBROUTINE ivreg(vsize,yin,xin,betaout,Vresout,E2out,out_err)
  USE Commonvars

  IMPLICIT NONE
  INCLUDE 'mpif.h'
  
  INTEGER, INTENT(IN)  :: vsize  
  REAL(8), INTENT(IN)  :: yin(vsize), xin(vsize,Nvarreg)
  REAL(8), INTENT(OUT) :: Vresout,E2out,betaout(Nvarreg),out_err
  REAL(8), ALLOCATABLE :: res(:)
  REAL(8)              :: Mres, xpx(Nvarreg,Nvarreg), xpy(Nvarreg), invxx(Nvarreg,Nvarreg)
  INTEGER              :: info, ii, jj
  CHARACTER            :: uplo='U'

  out_err = 0d0

  ALLOCATE(res(vsize))
  res = 0.0d0
  
  xpx = 0.0d0
  xpy = 0.0d0
  
  xpx = MATMUL(TRANSPOSE(xin(:,:)),xin(:,:))/DBLE(SIZE(yin))
  xpy = MATMUL(TRANSPOSE(xin(:,:)),yin(:))/DBLE(SIZE(yin))
  
  invxx=0.0d0
  DO jj = 1, Nvarreg
     invxx(jj,jj)=1.0d0
  END DO
  
  CALL DPOSV(uplo,Nvarreg,Nvarreg,xpx,Nvarreg,invxx,Nvarreg,info)
  
  IF (info /= 0) THEN
     DO jj = 1, Nvarreg
        IF (myrank == 0) WRITE(*,*) 'xpx', xpx(jj,:)
     ENDDO
     WRITE(*,*) 'MATRIX COULD NOT INVERT simulatemun', 'info', info, 'myrank', myrank, 'Nvarreg', Nvarreg
     DO ii = 1, vsize
        IF (myrank == 0) WRITE(*,*) 'y', yin(ii), 'x', xin(ii,:)
     END DO

     out_err = 1d0

  END IF
  
  betaout = MATMUL(invxx,xpy)

  res = yin(:) - MATMUL(xin(:,:),betaout(:))
  
!!$  Rsquared = (SUM((yin(:)-SUM(yin(:))/SIZE(yin(:)))**2d0) - SUM(res(:)**2d0)) / SUM((yin(:)-SUM(yin(:))/SIZE(yin(:)))**2d0)
!!$
!!$  if (Nvarreg==6) write(*,*) 'Rsquared', Rsquared

  Mres = SUM(res)/DBLE(vsize)
  Vresout = SQRT(SUM((res(:)-Mres)**2d0)/DBLE(vsize))
  E2out = SUM(res(:)**2d0)/DBLE(vsize)
  
  DEALLOCATE(res)
  
END SUBROUTINE ivreg
